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Abstract 

We present the results of molecular dynamics simulations of the extended simple point charge 
(SPC/E) model of water to investigate the thermodynamic and dynamic properties of stretched 
and supercooled water. We locate the liquid-gas spinodal, and confirm that the spinodal pressure 
increases monotonically with T, supporting thermodynamic scenarios for the phase behavior of 
supercooled water involving a "non-reentrant" spinodal. The dynamics at negative pressure show 
a minimum in the diffusion constant D when the density is decreased at constant temperature, 
complementary to the known maximum of D at higher pressures. We locate the loci of minima of D 
relative to the spinodal, showing that the locus is inside the thermodynamicaly metastable regions 
of the phase-diagram. These dynamical results reflect the initial enhancement and subsequent 
breakdown of the tetrahedral structure and of the hydrogen bond network as the density decreases. 
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I. INTRODUCTION 



Water is an important liquid in nature, and is also fundamental in chemical and technolog- 
ical applications. Although the individual water molecule has a simple chemical structure, 
water is considered a complex fluid because of its anomalous behavior |^, ^ |^, ^, |^. It 
expands on freezing and, at a pressure of 1 atm, the density has a maximum at 4°C. Addi- 
tionally, there is a minimum of the isothermal compressibility at 46°C and a minimum of the 
isobaric heat capacity at 35°C 0. These anomalies are linked with the microscopic structure 
of liquid water, which can be regarded as a transient gel — a highly associated liquid with 
strongly directional hydrogen bonds 0, ^. Each water molecule acts as both a donor and 
an acceptor of bonds, generating a structure that is locally ordered, similar to that of ice, 
but maintaining the long-range disorder typical of liquids. Despite the extensive work that 
has been done on water, many aspects of its behavior remain unexplained. 

Several scenarios have been proposed to account for the the anomalous behavior of the 
thermodynamic response functions on cooling, each predicting a different behavior for the 
liquid spinodal, the line of the limit of stability separating the region where liquid water is 
metastable from the region where the liquid is unstable, (i) According the stability-limit 



conjecture ||^, |T0[, the pressure of the spinodal line should decrease on cooling, become neg- 
ative, and increase again after passing through a minimum. It reenters the positive pressure 
region of the phase diagram at a very low temperature, thereby giving rise to a line of 
singularities in the positive pressure region, and consequently the increase in the thermo- 
dynamic response functions on cooling in the anomalous region is due to the proximity of 



this reentrant spinodal. (ii) The critical point hypothesis ||TT| , |12|, [l^, |T5|, |l^], proposes a 
new critical point at the terminus of a flrst-order phase transition line separating two liquid 
phases of different density. The anomalous increases of the response functions, compress- 
ibility, speciflc heat, and volume expansivity, is interpreted in terms of this critical point, 
(iii) The singularity-free hypothesis [§, 0, |l^ proposes that actually there is no divergence 
close to the anomalous region; the response functions grow by lowering temperature but re- 
main flnite, attaining maximum values. Both hypotheses are consistent with a non-reentrant 
spinodal, and so unambiguously identifying the spinodal can rule out at least on scenario. 

Water properties and anomalies can be strongly influenced by the physical or chemical 
properties of the medium ||, ^ |, [1^, The effect not only of applied pressure, but 
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also of negative pressure ("stretching") is remarkable. The study of the behavior of this 
fluid under negative pressures is relevant not only from the academic point of view, but 



also for realistic systems. For example, negative pressures are observed and seem to 
play an important role in the mechanism of water transport in plants. Therefore, properties 
that modify the structure of water, especially if this modification is similar to the effect of 
stretching (as is the case in some hydrogels |T^), also influence its dynamical behavior. 
Dynamic properties, such as the diffusion constant, have been studied in detail for water 



systems at atmospheric and at high positive pressures, both experimentally ^ as well 



as by computer simulations [^, |25|, ^ |29|, The increase of pressure increases 



the presence of defects and of interstitial water molecules in the network ||25|. They disrupt 
the tetrahedral local structure, weakening the hydrogen bonds, and thus increasing the 



diffusion constant However, a further increase in the pressure leads to steric effects 

which works in the direction of lowering the mobility. The interplay of these factors leads 
to a maximum in the diffusion constant ||2^, ^ at some high density pmax- Above this 
density (or corresponding pressure), the diffusion of water is in some sense like that of a 
normal liquid, controlled by hindrance, with the hydrogen bonds playing a secondary role. 
However, the behavior at very low p is less well understood. 

In this paper, we study how the thermodynamics and the dynamics of low-temperature 
water are affected by the decrease of the density. We perform molecular dynamics (MD) 
simulations of the SPC/E model of water in the range 210 K < T < 280 K and 
0.825 g/cm^ < p < 0.95 g/cm^. State points in this range have negative pressure, and 
are either a metastable stretched liquid, or a phase separated liquid-gas mixture. Additional 
details of the simulation procedure can be found in ref. pOf. Table | shows the thermody- 
namic and dynamic properties of the simulated state points as well as the time needed for 
equilibration and production runs. At the lowest T studied, the production runs are not 
long enough to reliably estimate D. 



II. LOCATION OF THE SPINODAL 



The effect of negative pressure is still not fully understood, and only a few experimental 

|33| , ^ |35[| have been performed on transport 



works |31|, |32[ and simulations |15|, |28| , |29|, |c 
properties of stretched water. In this negative pressure region of the phase diagram the sys- 
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tern is metastable, and becomes unstable beyond the spinodal line, so locating the spinodal 
we can ensure that our simulated state points lie in the metastable and not in the unstable 
region. Moreover, the shape of the spinodal can test the stability-limit conjecture against 
the critical point hypothesis and the singularity-free interpretation, so we first locate the 
density and pressure of the spinodal, which we denote Psp(T) and Psp(T), respectively. 

We consider the the P-p behavior along isotherms. For each T, there is a minimum in P 
at p ^ 0.86 g/cm^. At lower p, P increases, signaling cavitation of the liquid. For densities 
below psp, Kt = p^^{dp/dP)T becomes negative, so the liquid is unstable with respect to the 
vapor phase. Therefore, the locus of infinite compressibility locates the spinodal line, and 
sets the density below which our simulated state points are no longer thermodynamically 
stable. 

We confirm which state points are "single phase" and which are "phase separated" by 
calculating the oxygen-oxygen partial structure factor 

1 ^ 

Siq) = ttT.^-"''-^'''''''^ (1) 

Figure | shows S{q) at the the lowest T (T = 210 K) for each density simulated. At the 
lowest density, p = 0.85 g/cm^, S{q) for g ^ 10 nm~^ increases significantly, signaling the 
presence of large-scale structure in the system due to phase separation of liquid and gas. 
Similar behavior is observed for other state points that would appear to be in the unstable 
regime from Fig. [l| ||36|| . 

To test whether the minima in the P-p isotherms arise from cavitation of the liquid, we 
fit at each T a 5th-order polynomial to the P-p isotherms, and estimate Psp{T) and Psp{T) 
by the minimum of the polynomial fit. We find that psp{T) occurs in the narrow range 
0.853 < Psp < 0.874, and so the spinodal density is nearly T independent. Our results 
indicate that Psp(T) decreases monotonically with decreasing T (Fig. |^), ruling out the 
possibility of a re-entrant spinodal in the temperature range studied. 

We also compare our estimate of the spinodal with estimates from two other sources for 
the SPC/E potential: (i) an estimate of the upper bound on -Psp (7") and (ii) the spinodal 
estimated by an approximate analytic equation of state (EOS) [|16]. Figure | shows that 
our estimate of -Psp (7") is less than (and hence consistent with) the upper bound estimate 



of Ref. [27|. The spinodal obtained from the EOS in Ref. |T6| is quantitatively reasonable 



for T ^ 280 K, but decreases far more quickly with decreasing T (reaching ^ —2 GPa at 



T = 210 K) than our simulations indicate. This dramatic decrease hkely results from the 
fact that the fitting procedure requires P along a reference isotherm, which was performed 
at T = 300 K in the density range 0.95 < p < 1.4 g/cm^. This fit is not reliable outside of 
this density range, where the spinodal lies, causing an erroneous spinodal estimate. These 
errors are compounded as the fit is pushed further from the reference T. To improve the 
estimate, we repeat the fitting procedure using a new reference isotherm T = 210 K, and 
expand the density range down to p = 0.85 g/cm^. The results of this fit are shown in Fig. 
and we obtain more reasonable agreement with our direct simulation results. Ref. [|l^] also 
used this EOS to predict a liquid-liquid phase transition for T ^ 130 K; this prediction 
should not be affected by using data at lower p, since the transition lies within the density 
range of the original fits. 

III. DYNAMIC PROPERTIES 

Having established the limit of stability at negative pressures, we next analyze the dy- 
namics of stretched water. The effect of extreme conditions on the flow of the liquid we 
assess by calculating the diffusion constant D, defined by the asymptotic value, obtained 
by linear regression, of the slope of the mean square displacement versus time. We show D 
along isotherms in Fig. ^ For T < 260 K, D has a minimum value at p ~ 0.9 g/cm^, which 
becomes more pronounced at lower T (Fig. ^). This behavior can be understood consid- 
ering the structural changes that occur with decreasing density. At low T, the decreased 
density enhances the local tetrahedral ordering, which leads to a decrease in D. Further 
decreases in density reduces the stability of the tetrahedral structure and causes an increase 
ofD. 

The location of the minimum is near the ice Ih density ~ 0.915 g/cm^, which is the 
density where the perfect tetrahedral order occurs. The behavior of the minimum of D, 
D„un(T), complements the known behavior of -Dmax(^) for the same model [^, |3y, |35| , 
where a maximum occurs due to breaking hydrogen bonds at high pressure; the density of 
the Dmin(^) increases slightly with increasing T, while the density of D^.^^{T) decreases with 
increasing T This is expected, since the range of densities where anomalous behavior 
occurs expands with decreasing T. We show the loci of D^^i^iT), D^i^^iT), along with the 
spinodal and the temperatures of maxima density, Tmd{P) in Fig. |. 
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IV. CONCLUSIONS 



Water exhibits a very complex structure and its properties and anomalies are strongly 
influenced by variations of pressure. For high densities (p > Pmax), water behaves as a 
normal liquid and the decrease of D with increasing pressure is governed by steric effects. For 
Pmin < P < Pmax, as the pressure is decreased, the presence of defects and interstitial water 
decrease, the tetrahedral structure dominates, with stronger hydrogen bonds. This process 
reaches its maximum at p = pmin ~ Pice- Further stretching destabilizes the hydrogen bond 
network, leading to an increase in mobility. The locus of Dmin roughly tracks the spinodal, 
not surprising since the same breakdown of tetrahedral order that gives rise to Dmin also 
facilitates cavitation. 
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TABLE I: Temperature, density, diffusion constant, potential energy, pressure, equilibration t, 
and production tprtiicnes. The uncertainty in D is roughly ±0.008. 
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FIG. 1: Dependence of the pressure on the density, along isotherms, for the simulations reported 
here. 
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FIG. 2: At r = 210 K, the structure factor S{q). Note the increase of S{q) at small q aX p = 
0.85 g/cm^, indicating phase separation. 
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FIG. 3: Location of the spinodal (o), showing no evidence for re-entrant behavior. For comparison, 
we also show the upper bound on the spinodal (filled □) from ref. [^] and the estimated spinodal 
(dotted line) from the analytic free energy expression of ref. . 
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FIG. 4: (a) Dependence of the diffusion constant D on p along isotherms (for p < 1.0 g/cm'^). 



Open symbols are the new simulations we report, and filled symbols are from ref. |2£]. The dotted 
line separates liquid state points from phase separated state points, but is not an indication of the 
exact psp{T), which varies slightly with T. (b) Full p dependence of D, also showing the maxima. 
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FIG. 5: Relation of the loci of maxima and minima of D with Tmd{P) and the spinodal. Open 
symbols are from the present work, and filled symbols are from ref. [^]. There is no maximum 
and minimum of D for temperatures T > 280 K Fig. ^ 

and ref. [|| 
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